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I.       INTRODUCTION 

A.  DESCRIPTION  OF  THE  PROBLEM 

The  study  of  dynamic  responses  in  a  submersible  vehicle  using  a  nonlinear 
analysis  is  important  in  determining  operating  envelopes  for  the  vehicle.  Previous  studies 
have  shown  that  in  straight  line  motion,  the  coupling  of  the  sway,  yaw,  and  roll  equations 
produce  oscillating  losses  of  stability  in  the  system.  [Ref.  5]  Introducing  a  nonzero  pitch 
angle  to  the  vehicles  motion  will  allow  us  to  study  the  changes  to  the  stability  domain  for 
a  variety  of  environmental  conditions. 

In  our  work,  the  linearized  equations  of  motion  are  studied  using  an  eigenvalue 
analysis  to  determine  the  systems  stability  through  a  variety  of  operational  parameters.  A 
nonlinear  analysis  is  then  conducted  to  assess  the  stability  characteristics  of  the  resulting 
limit  cycles  and  their  impact  on  the  operating  characteristics  of  the  vehicle. 

B.  OBJECTIVES  AND  OUTLINE 

In  this  thesis  we  expand  on  previous  thesis  work  which  examined  the  problem  of 
stability  in  straight  line  motions  of  a  submersible  vehicle  [Ref.  13].  The  primary  cause 
for  this  loss  of  stability  is  the  coupling  between  the  sway/yaw/roll  equations  of  motion  for 
a  submersible  vehicle.  We  know  that  the  loss  of  stability  creates  stable  limit  cycles  in 
straight  line  motion.  This  work  analyzes  the  effects  of  introducing  a  nonzero  pitch  angle 
to  the  generic  equations  of  motion  in  order  to  determine  their  effect  on  the  creation  of 
limit  cycles. 


The  model  used  for  this  work  is  a  variant  of  the  Swimmer  Delivery  Model  used  in 
[Ref.  5]  for  which  a  generic  set  of  hydrodynamic  and  geometric  properties  are  available. 


II.        EQUATIONS  OF  MOTION 

A.  COORDINATE  SYSTEM 

A  moving  coordinate  system  was  used  for  our  analysis  with  the  origin  fixed  at  the 
vehicles  center  of  buoyancy.  The  x-axis  is  fixed  to  the  longitudinal  plane  of  symmetry  for 
the  vehicle,  the  y-axis  is  positive  to  starboard,  and  the  z-axis  is  positive  downward.  All 
symbols  used  in  the  development  of  the  equations  of  motion  are  summarized  in  Table  1 . 

B.  GENERAL  FORM  OF  THE  EQUATIONS  OF  MOTION 

The  equations  of  motion  for  a  submersible  vehicle  in  the  horizontal  plane  are 
written  as  follows: 

Sway  equation: 

m[v  +  Ur-wp  +  xc(pq  +  r)-yG(p2  +r2 )+  zG (qr -  p)]= 
Y.p  +  Y,r  +  Ypqpq  +  Yqrqr  +  YvUv  +  Ymvw  + 

Y5r  U  28r  +  F,v  +  YpUp  +  YrUr  +  Yvqvq  +  Ywp  wp  +  Ywr  wr  + 
(W-5)cos0sin0- 

f""  [CD  h(xXv  +  xrf  +  CD  b(x\w  -  xqf  \V     *\  dx 

Jx«»'        y  -  u cf{x)  U; 


Yaw  equation: 


IJ  +  {lZz-Iyy)p<l-IXy(p2-42)-Iyz(Pr  +  4)+Ixz(4r-p)+ 

m[xG  (v  +  Ur-  wp)-  yG  if)  —  vr  +  wq)\  = 

Npp  +  N,r  +  Npqpq  +  Nqrqr  +  Nl,Uv  +  Nmyw  + 

NsU28r+Nvv  +  NpUp  +  NrUr  +  Nvqvq  +  Nwpwp  +  Nwrwr  + 
(xcW-xBB)cosdsm(P  +  (ycW-yBB)sme  +  U2Nprop- 

P"  [cD  h{xlv  +  xrf  +  CD  b{x\w  -  xqf  fc^i xdx 

Jx....i  y  J 1     I  v  1 


uM) 


(2) 


Roll  equation: 

m\yG (w-Uq  +  vp)-  zG {v  +  Ur-  wp)]  = 

Kpp  +  K,r  +  Kpqpq  +  Kqrqr  +  KvUv  +  Kvwvw  + 

KnropU2  +  Kvv  +  K  Up  +  KrUr  +  Kvq  +  Kwp  +  Kwrwr  + 


prop 


(yGW  -  yBB)cos6 cos0  -  (zGW  -  zBB)cos6 sin0 


(3) 


The  rotational  velocity  equation  around  the  x-axis: 


<j)  =  p 


(4) 


UCf  denotes  the  cross-flow  velocity: 


Ucf  (x)  =  yj(v  +  xrf  +  (w  -  xq)2 


(5) 


Cn 

quadratic  drag  coefficient 

dr 

rudder  deflection 

h(x) 

local  height  of  hull 

( *■  xx  >  *  yy  *  zz  ' 

vehicle  mass  moments  of  inertia  about  body  axes 

(■*  xy>*  yz>*  zx  ) 

cross  products  of  inertia 

(K,M,N) 

moment  components  along  the  three  axes 

m 

vehicle  mass 

(p,q,r) 

rotational  velocity  components  along  the  body  axes 

(f,q,y) 

Euler  angles 

U 

constant  vehicle  speed  along  the  x  -axis 

(u,v,w) 

translational  velocities  about  (x,y,z)  axes 

(x,y,z) 

distances  along  the  three  body  axes 

(X,Y,Z) 

force  components  along  the  body  axes 

(xG,yG,zG) 

coordinates  of  the  center  of  gravity 

(xB,yB>zB) 

coordinates  of  the  center  of  buoyancy 

■^  nose 

fore  coordinate  of  vehicle  body 

x  tail 

aft  coordinate  of  vehicle  body 

Table  1:  Nomenclature 

C.       SIMPLIFICATIONS 

We  must  simplify  the  equations  of  motion  in  order  to  reflect  the  fact  that  we  are 
analyzing  motion  about  the  y-axis.  The  simplifications  that  we  employ  are: 


Acceleration,  w ,  in  the  z-direction  is  zero. 


Acceleration  in  the  longitudinal  direction,  u  ,  is  zero. 

Rotational  velocity,  q,  and  acceleration,  q ,  in  the  y-  direction  are  zero 


D.       SIMPLIFIED  EQUATIONS  OF  MOTION 

After  applying  the  above  simplifications,  the  equations  of  motion  become: 

Sway  equation: 

m[v  +  Ur-wp  +  xcr-yc(p2+r2)+zGp]= 

Ypp  +  Yrr  +  YvUv  +  Yvwvw  +  Y5U2dr+Yvv  +  YpUp  +  YrUr  +  Ywpwp  +  Ywrwr  + 

(W-£)cos0sin0- 

f"  [cDh{xlv  +  xr)2  +CDb{xW}^\dx  (6) 

Jx""'  "  Ucf\X) 


Yaw  equation: 


^r-lxyP2-Iyzpr  +  Ixzp  + 

m\xc  (v  +  Ur-  wp)-  yG  if)  -  vr)\= 

N.p  +  N,r  +  NvUv  +  Nvwvw+N5U2Sr+Niv  +  NpUp  +  NrUr  +  Nwpwp  + 

Nwrwr  +  {xcW-xBB)cosdsm(t)  +  {ycW-yBB)smd  +  U2Nprop- 

£j  [cDh(x\v  +  xrY  +CDMxy}^xdx 


M?  (7) 


Roll  equation: 


txxP  +  txyPr-^y-txzr  +  mhGM-Zciv  +  Ur-wp)^ 

KpP  +  Krr'+KvUv  +  Kvwvw+KpropU2+Kvv+KpUp  +  KrUr  + 

K    wp  +  Kwrwr  +  (yGW  -  y  BB)  cos  6  cos  (f)-(zGW  -  zBB)cos6s'm  (j)  ^ 


Roll  rate: 


<p  =  p 


(9) 


III.       LINEAR  ANALYSIS 


A.       LINEARIZATION 


The  simplified  equations  of  motion  can  be  written  in  the  matrix  form: 


Ax  =  Bx  +  g(x) 


(10) 


where  the  state  vector,  x  ,  is  defined  as, 


and  the  state  matrices  are, 


x  = 


A  = 


m-Y, 

mxG-Y. 

-mzc-YP 

0 

mxG-Ni 

/«-*, 

-I~-N* 

XZ                   P 

0 

-mzG-Ki 

Jxz  ~Kr 

XX                p 

0 

0 

0 

0 

1 

and 


B  = 


Yv.  +  Yww 

KvU  +  Kww 
0 


YrU+Ywr-mU 

-mxGU  +  NrU  +  Nwrw 

mzGU  +  KrU  +  Kwrw 

0 


YpU  +  Ywpw  +  mw  0 


NpU+Nwpw 

0 

mzGw+KpU  +  Kwpw 

0 

1 

0 

The  g(x)  term  remains  as  given  [Ref.  13]. 

The  nonlinear  terms  are  linearized  around  a  nominal  point: 

xo  =[vo>ro>Po><l)o]T  =° 

A  Taylor  series  expansion  is  applied  to  the  nonlinear  terms  about  the  nominal  point,  Xo  , 
and  keeping  only  the  linear  components,  the  equations  of  motion,  written  in  matrix  form, 
become: 


A'x  =  Bx 


(11) 


where, 


and 


A'  =  A 


B'  = 


YVU 

YrU-mU 

YS> 

(W-B)cosd 

NVU 

-mxGU  +  NrU 

NpU 

(xGW  -xBB)cos6 

KVU 

mzGU  +  KrU 

KPu 

{-zgW  +  zbB)cos6 

0 

0 

1 

0 

To  assess  the  dynamic  stability  of  the  vehicle,  an  Eigenvalue  analysis  is  performed  in  the 


next  section. 
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B.       LOSS  OF  STABILITY 

An  Eigenvalue  analysis  is  used  to  determine  the  stability  of  the  linearized  system. 
Stability  is  dependent  on  the  location  of  the  Eigenvalues  or  the  roots  of  the  characteristic 
equation: 

det(£'-yU'  =  0)  (12) 

in  the  polynomial  form: 

AA4+5A3+CA2+Z)A  +  £  =  0  (13) 

The  coefficients  equation  (13)  are  given  in  [Ref.  13,  pg.  11].  Using  Routh's  criterion  we 
can  examine  the  stability  of  the  system.  The  following  conditions  must  be  applied  to  the 
characteristic  equation  (13)  to  ensure  all  roots  have  negative  real  parts: 

BCD-AD2  -EB1  >0  (14) 

£>0  (15) 

If  E  is  less  than  zero,  one  real  root  of  ( 13)  becomes  positive  and  the  system  will  become 
unstable  in  a  divergent  manner  [Ref.  9].  This  is  the  case  of  a  directionally  unstable  ship 
which  is  well  known  in  the  literature  [Ref.  3].  If  the  condition  in  (14)  is  not  met,  it  means 
that  there  is  at  least  one  complex  conjugate  root  with  real  parts  and  will  result  in  an 
oscillatory  response,  indicating  loss  of  stability. 
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To  analyze  the  limiting  case  of  loss  of  stability  equation  (14)  must  be  solved  in  the 
following  form: 

BCD- AD2  -EB2  =0  (16) 

The  result  of  this  equation  will  be  a  curve  of  zg  as  a  function  of  xg  and  will  be  our  locus 
for  loss  of  stability.  We  can  express  the  coefficients  of  equation  (16)  in  the  form: 

A  =  AlzG+A2zG+A3  (17) 

B  =  Blzl+B2zG+B3  (18) 

A\,  A2,  and  A3  are  of  the  form  given  in  [Ref.  13,  pg.  14]. 

B\  and  B3  are  of  the  form  given  in  [Ref.  13,  pg.  14].  With  the  addition  of  a  pitch  angle,  w, 


B2  takes  the  form: 

B2  =-m{KvU\lx  -N>)-m{la  -NrppU) 
+  mY.  (UNr  - UmxG  )+  mK,  (UNr  - Umxc  ) 
-m{Np  +lJ{JYr  -Um)+m{N pU\mxG  -Yr) 
~m{Kr  -Ixz  XNvU)+mUYp{mxG  -Nv) 
-mU{Np  +Ixz\m-Yv)+mUKr(mxc  -Nv) 
+  mzG  w(m  -YiXlzz-Nr )-  mzG  w(mxG  -  Y,  \mxG  -  Nv ) 

C  =  C}zl+C2zG+C,  (19) 
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C\  and  Cz  are  of  the  form  given  in  [Ref.  13,  pg.  15].  With  the  addition  of  a  pitch  angle, 
w,  C2  takes  the  form: 

C2  =  mU(mxG  -NvlYpU)-mUYp(NvU)-mUKr(NvU) 
+mU{Np+Ixz\YvU)-mU{NpuXm-Yi) 

+  W{m-Yv\la  -Nr)-W{mxG  -YrtmxG  -N,) 
-m(XBB-xGW\mxG  -Yr)+m{UNr  -UmxG^YpU) 
-m(UYr  -Um^N pU)+m{KvU\UNr  -UmxG) 
-mz0w{m-Yv\UNr  -UmxG)-mzGw(YvUlla  -Nr) 
-mzGw(mxG  -Y-\NvU)+mzGw{mxG  -N-\UYr  -Um) 

D  =  DxzG+D2  (20) 

D\  is  of  the  form  given  in  [Ref.  13,  pg.  16].  With  the  addition  of  a  pitch  angle,  w,  D2 
takes  the  form: 

D2=UKr{xBB-xGWlm-Y,)-UKr(NvUiYpU) 
+UKr{NpU\YvU)-(Kf-IxzXxBB-xGwiYpU) 

+  Ki{xBB-xGWlUYr-Um)-(KvU\xBB-xGWXmxG-Yr) 
+  {KvUXUNr  -UmxGiYpu)-{KvU\UYr  -Un^N plj) 
-{KpU\YvUlUNr-VmxG)+{KpU\UKr-UmlNvU) 
+  mzG  w(YvU  \UNr  - UmxG  )-  mzg  w(UKr  -  Um\N \,U )  ' 

The  equation  for  the  coefficient,  E,  remains  unchanged  [Ref.  13,  pg.  16].  Applying  the 
stability  criterion,  equation  (16),  and  utilizing  them  in  the  resulting  fifth  order 
polynomial,  F,  [Ref.  13,  pg.  18]  we  are  able  to  solve  F  using  the  MATLAB  program  in 
Appendix  A.  The  curves  show  zg  as  a  function  of  xq  and  we  show  results  for  varying 
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pitch  angles,  w,  and  varying  forward  velocities,  U.  On  all  of  the  following  graphs  the 
pitch  angle  is  varied  from  10  degrees  to  -10  degrees  in  increments  of  five  degrees.  The 
top  curve  is  the  10  degree  curve  and  the  bottom  curve  is  the  -10  degree  curve. 
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Figure  1:  ZG  vs.  XG  for  U  =  2  ft/sec 
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Figure  2:  ZG  vs.  XG  for  U  =  3.5  ft/sec 
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Figure  3:  ZG  vs.  XG  for  U  =  5  ft/sec 
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Figure  4:  ZG  vs.  XG  for  U  =  6.5  ft/sec 
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Figure  5:  ZG  vs.  XG  for  U  =  8  ft/sec 
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From  the  results  of  Figures  1  through  5,  the  following  conclusions  can  be  drawn: 

1.  In  all  cases,  a  sufficiently  high  metacentric  height  is  required  in  order  to  ensure 
vehicle        stability.  Regions  of  parameters  that  fall  below  the  critical  curves 
correspond  to  dynamic  instability. 

2.  The  critical  metacentric  height  that  is  required  for  dynamic  stability  is  an  increasing 
function  of  both  vehicle  speed  and  trim  angle. 

3.  Static  stability  alone  does  not  necessarily  ensure  dynamic  stability  of  motion  during 
the  turn. 

4.  The  loss  of  stability  experienced  here  is  a  dynamic  loss  of  stability.  At  the  critical 
metacentric  height,  one  pair  of  complex  conjugate  eigenvalues  possesses  a  zero  real 
part.  This  is  an  oscillatory  loss  of  stability,  which  can  not  be  predicted  by  considering 
the  uncoupled  sway/yaw  equations  of  motion  alone. 
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IV.       NONLINEAR  ANALYSIS 

A.  INTRODUCTION 

In  the  previous  chapter  we  performed  a  linear  analysis  of  the  equations  of  motion, 
observing  that  with  changes  to  specific  parameters  (jcg,  Zg,  w)  it  is  possible  to  pass 
through  a  stable  region  to  a  region  of  loss  of  stability. 

In  previous  work  [Ref  13]  it  was  shown  that  these  bifurcations  to  periodic 
solutions  were  all  supercritical.  This  means  that  limit  cycles  were  produced  after  the  loss 
of  stability.  By  introducing  a  pitch  angle,  w,  in  the  nonlinear  analysis  we  will  analyze 
whether  the  bifurcations  to  periodic  solutions  will  remain  supercritical  and  what  changes 
occur  to  the  limit  cycles  themselves. 

B.  THIRD  ORDER  EXPANSIONS 

Our  linearized  system  was  written  in  the  form  of  equation  (11),  ignoring  the 
nonlinear  terms.  Including  the  nonlinear  terms  changes  the  form  of  equation  (1 1)  to:  . 

A'x  =  B'x  +  g(x)  <21> 


where: 
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gW= 


82(x) 

sAx) 
sAx) 


Keeping  terms  up  to  third  order,  the  matrix  g(x)  can  be  written  in  the  vector  form: 


g(x)=gV(x)+gV(x)+c(x) 


(22) 


where  g^2\x)  contains  the  second  order  nonlinear  terms  and  g^{x)  contains  the  third 
order  nonlinear  terms.  The  cross  flow  integrals  and  the  second  order  nonlinear  terms 
remain  unchanged  with  the  addition  of  a  pitch  angle,  w.  [Ref.  13,  pg.22,  23]  However 
the  third  order  nonlinear  terms  take  the  form: 


where: 


g«= 


g(3)=_7(3) — (w  -  By  cosd 
6 

g(3)=_/(3)_I( xW-xBBy  cosd 
o 

g?=UzgW-zBBy  cose 
o 


The  Taylor  series  expansion  yielding  the  second  and  third  order  linear  terms  of  the  cross 
flow  integrals,  and  the  inverse  of  the  system  matrix,  (A')~  ,  remain  unchanged.  [Ref.  13, 
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pg.  24,  25].  However  the  B'  matrix  has  changed  as  shown  in  the  Linear  Analysis  section 
of  this  work.  These  changes  to  the  4  x  4  F  matrix  are: 


F  = 


Ml  M2 


D 
0 


D 
0 


D 

1 
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£> 

£> 

D 

D 

M>. 

F 

M* 

M* 

D 

D 

D 

D 

^L 

F* 

*k 

5i 

0 


where: 


F11=fl11(7vJ7)+a12(iVvf/)+a13(^vt7) 

F12  =au(Yr  -mp+a12(Nr  -mxG)U  +  an(Kr  +  mzG)U 

Fn=au(YpU)+an{NpvU)+an{KpvU) 

M4  =  an  (w  _  #)cos0  +  a12  (*GW  -  Jts£)cos0  +  a13  (-  zGW  +  zBB)cosd 

F*  =a2l{YvU)+a22{NvU)+a23{KvU) 

M22   =  a2i  (m  -^X7  +  «22  (#v  _m*G  X7  +fl23  (^v  +  m^G  X7 

m*  =a2l{YpU)+a22{NpvU)+a23(KpvU) 

F24  =  a21  (W  -  B)cosO  +  a22  (xcW  -xBB)cos6  +  a2l  (-  zcW  +  zBB)cosd 

Fn  =a3l(YvU)+a,2{NvU)+a„(KvU) 

F32  =  a31  \YT  -  m)U  +  <a32  (Nr  -  mxG  )(J  +  a33  (Kr  +  mzG  ]U 

m*  =a31{Ypu)+a,2{NpU)+a33{KpU) 

FM  =  a31  (W  -  Z?)cos  0  +  a32  (xG  W  -  xBB)cos  6  +  a33  (-  zG  W  +  zBB)cos  6 
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The  remaining  elements  of  the  nonlinear  and  constant  terms  remain  unchanged  [Ref .  13, 
pg.  28,  29]. 

C.  COORDINATE  TRANSFORMATIONS 

To  continue  our  analysis  it  is  necessary  to  bring  our  transform  our  coordinate 
system  from  state  space  to  a  normal  coordinate  system.  This  transformation  is  performed 
in  the  manner  given  [Ref.  13,  pg.  29,  30]. 

D.  CENTER  MANIFOLD  EXPANSIONS 

The  center  manifold  expressions  are  of  the  form  given  [Ref.  13,  pg.  30-33]. 
However,  the  coefficients  %/=!,  2,  3;j=5,  6,  7  are: 


,         au        (    2         2  \ 


D 


+  -j^VXyml\  +Iyz^xm2X+yGm2ymu) 


D 


rIjym3lm2l  +Iyzmlx  +myGmum3l  N 
+  y„Wm24,-yBBm24] 


(23) 


h,6  =~WyG  (2m31  m32  +  2m2\  +  W22  ) 


a 


+  ■ 


12 


D 


"13 

D 


2I39m3lm32  +I,Z  (mMm22  +  m32m2l ) 
+  yG(m2]mn+m22mn) 


I  x\  (m3im22  +^32m21  )+^  yz171 


yz'"'2\m22 


+  myG  {mumn  +  ml2m3l )+  2{ycW  -  yBB)mAXm42 


(24) 
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j  ^11  (l  2   \ 

Kn  =~yc\m32  +m22) 


D 


+  -^[1^32  +Iyzm32m22+ycm22mn) 

-^■[IXym32m22  +/vzm22  +^GW12W32  +  ^  GW  "  ^5K2  J 
1  ^21  (l  2   \ 

h,5=  —  yG^+m2x) 

+  -j^{lXym3l  +IyzmMm2l+yGm2^n) 


a 


23 

D 


flXym3Xm2\  +Iyzml\  +myGmiim3l  ^ 


+  yGWml1-yBBmli 


(25) 


(26) 


h,6  =-^-yG(2m3im32+2m2l+m22) 


+  ■ 


D 


a 


23 


D 


2Ixym3im32  +/y2(m3im22  +  ™32W21  ) 

+  yG(m2lml2+m22mu) 

Ix>(mnm22  +mnmlx  )+2Iyzm21m22 

+  myG(mum32+mnm2i)+2(ycW-yBB}n4]m42 


(27) 


1  ®2\  (2  2   \ 

l2,l   =  —  yc{m32+m22) 


D 


+  -T-(/^m322  +Iyzm32m22+yGm22mn) 


D  yxy 

^-l7xvm32m22  +Iyz^22  +  "V Gm\2m32  +  (^GW  "  ^KL  J 

7  ^31  f      2  2   ^ 

+  ^r{lXym32l  +Iyzm3im21+yGm21*nn) 


(28) 
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— ^r(7^m3im2i  +IyZml\  +rnyGmnmll+yGWml,  -yBBml\) 
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a. 


h,6   =~ftyG  (2m31  m32  +  2m21  +  W22  ) 


+  ■ 


'32 


D 
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D 


2Ixym3lmn  +  Iyz  (m3lm22  +m32m2i ) 
+  yc(m2lml2+m22mu) 

/^(m3im22  +mnm2x)+2Iyzm2Xm22 

+  myG(mum32+mnmil)+2(yGW-yBB)m4lm 


a3\  (      2  2    \ 

h,l   =  —  yG\m32+m22) 

+  -^-{lxymV  +Iyzmnm22+yGm22mn) 


D  L 


xsmi2m22 


I   ml2+myGml2m32  +{yGW- yBB)m242] 


(29) 


(30) 


(31) 


E.       AVERAGING 

The  procedure  for  averaging  the  equations  up  to  the  third  order  is  the  same  one 
given  in  [Ref.  13,  pg.  37-40].  The  addition  of  a  nonzero  pitch  angle  yields  new 
coefficients  ly,  I=\,  2,  3; 7=1,  2,  3,  4.  They  are: 
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3E0mumn  +  3Ex[mx22m2x  +  2ranra12ra22)+3£'2(ra22mII  +  2m2xm22mX2) 
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The  remaining  procedure  for  determining  the  equation  for  the  radius  of  the  resulting  limit 
cycles  is  identical  to  the  one  described  in  [Ref.  13,  pg.  43,  44],  which  is: 


R  =  a'eR  +  KR: 


(43) 
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F.        LIMIT  CYCLE  ANALYSIS 

At  steady  state,  R  =  0 ,  equation  (43)  becomes: 

0  =  R(a'e  +  KR2)  (44) 


Equation  (66)  has  two  solutions: 


R  =  0  (45) 


*-tT  (46) 


Equation  (45)  is  the  trivial  solution  and  gives  no  useful  information.  Equation  (46)  gives 
a  constant  amplitude  limit  cycle,  R,  in  the  cartesian  coordinate  system.  This  limit  cycle 
will  exist  if  the  quotient  inside  the  radical  sign  is  positive: 


—  a  e 
K 


>0  (47) 


This  condition  is  necessary  for  R  to  be  a  real  number.  Since  a '  is  always  negative,  the 
existence  of  limit  cycles  depends  on  the  parameter  K.  The  introduction  of  a  nonzero  pitch 
angle  does  not  change  the  dependence  the  limit  cycle  has  on  the  parameter  K.  Stated,  this 
dependence  is: 

•  If  K<0,  periodic  solutions  exist  and  they  are  stable. 

•  If  K>0,  periodic  solutions  exist  and  they  are  unstable. 
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G.       RESULTS  AND  DISCUSSION 

Results  for  the  stability  parameter,  K,  are  presented  in  Figures  6  through  15 
Figures  6  through  10  provide  results  for  K  at  a  constant  forward  speed,  U,  and  varying 
pitch  angles.  The  pitch  angles  were  varied  from  positive  10  degrees  to  negative  10 
degrees  in  5  degree  increments.  In  Figures  6  through  10,  the  bottom  curve  represents 
solutions  for  positive  10  degrees  and  the  top  curve  represents  negative  10  degrees.  It  is 
clear  that  all  values  of  K  are  negative  indicating  they  are  stable  solutions.  Notice  that 
decreasing  pitch  angles  the  solutions  for  K  become  less  negative,  indicating  that  while 
they  are  stable,  these  solutions  are  less  stable  than  those  at  the  higher  pitch  angles.  In 
Figures  1 1  through  15,  solutions  for  the  stability  parameter,  K,  are  again  represented,  but 
with  the  pitch  angle  held  constant  and  the  forward  speed,  U,  varied  from  2  ft/sec  to  8 
ft/sec  in  1.5  ft/sec  increments.  The  bottom  curve  in  Figures  1 1  through  15  represents  U  = 
2ft/sec.  As  forward  speed  increases,  the  curves  representing  the  stability  parameter  K 
tend  to  become  more  negative,  but  in  a  more  pronounced  fashion  than  those  where  U  was 
held  constant. 
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Figure  6:  XG  vs.  K*GAMMA  for  U=2  ft/sec 
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Figure  7:  XG  vs.  K*GAMMA  for  U=3.5  ft/sec 
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Figure  8:  XG  vs.  K*GAMMA  for  U=5  ft/sec 
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Figure  9:  XG  vs.  K*GAMMA  for  U=6.5  ft/sec 
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Figure  10:  XG  vs.  K*GAMMA  for  U=8  ft/sec 
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Figure  11:  XG  vs.  K*GAMMA  for  THETA  =10  deg 
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Figure  12:  XG  vs.  K*GAMMA  for  THETA=5  deg 
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Figure  13:  XG  vs.  K*GAMMA  for  THETA=0  deg 
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Figure  14:  XG  vs.  K*GAMMA  for  THETA=-5  deg 
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Figure  15:  XG  vs.  K*GAMMA  for  THETA=-10  deg 
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V.       CONCLUSIONS  AND  RECOMMENDATIONS 

This  thesis  presented  a  continuing  study  of  the  formation  of  limit  cycles  due  to  the 
coupling  of  the  sway/yaw/roll  equations  of  motion.  We  have  shown  that  loss  of  stability 
occurs  in  the  form  of  stable  limit  cycles  and  that  the  addition  of  a  nonzero  pitch  angle  does 
not  significantly  affect  the  formation  of  these  limit  cycles.  Through  a  linear  analysis  of  the 
sway/yaw/roll  equations  of  motion  we  demonstrated  that  the  addition  of  a  nonzero  pitch 
angle  affects  the  domain  of  stability  of  straight  line  motion,  especially  at  higher  speeds. 
This  was  validated  by  the  nonlinear  analysis  as  well.  As  a  recommendation  for  further 
study  in  this  area  we  suggest  that  the  analysis  be  expanded  to  include  coupling  into  the 
heave/pitch  directions  of  motion  as  well  as  the  effects  automatic  path  control. 
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APPENDIX 

The  following  is  a  list  of  the  MATLAB  and  FORTRAN  programs  used  in  this  thesis. 
Complete  printouts  of  the  programs  accompany  this  list. 

•  THESIS1.M 

A  MATLAB  program  for  performing  the  linear  analysis  section  of  this  thesis. 

•  HOPF_lNEW.FOR 

A  FORTRAN  program  for  performing  the  nonlinear  analysis  section  of  this  thesis. 
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%  THESIS l.M 

% 

%  LOSS  OF  STABILITY 

clear 

a=l; 

W=  12000; 

DCX=1760;  IYY=9450; 

IZZ=  10700;  DCZ=0;  DCY=0; 

IYZ=0;  L=17.425;  RHO=1.94; 

G=32.2;  U=8.0;  M=W/G;  B=W; 

THETA=0; 

THETA=THETA*pi/l  80; 

OMEGA=U*tan(THETA); 

ND1=0.5*RHO*LA2; 

%  DEFINE  HYDRODYNAMIC  COEFFICIENTS 


YPDOT=  1 .270e-04*ND  1  *LA2; 

YVDOT=-5.550e-02*NDl*L; 

YRDOT=  1 .240e-03*ND  1  *LA2; 

YP=3.055e-03*NDl*L; 

YPOMEGA=0; 

YP=YP*U+YPOMEGA*OMEGA+M*OMEGA; 

YV=-9.310e-02*NDl; 

YVOMEGA=0; 

YV=YV+YVOMEGA*OMEGA; 

YR=-5.940e-02*NDl*L; 

YROMEGA=0; 

YR=YR+YROMEGA*OMEGA; 

NPDOT=-3.370e-05*NDl*LA3; 
NVDOT=  1 .240e-03*ND  1  *LA2; 
NRDOT=-3 .400e-03  *ND  1  *LA3 ; 
NP=-8.405e-04*NDl*LA2; 
NPOMEGA=0; 

NP=NP+NPOMEGA*  OMEGA; 
NV=-1.484e-02*NDl*L; 
NVOMEGA=0; 

NV=NV+NVOMEGA*OMEGA; 
NR=- 1 .640e-02*ND  1  *LA2; 
NROMEGA=0; 
NR=NR+NROMEGA*OMEGA; 
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KPDOT=-1.01e-03*NDl*LA3; 

KVDOT=  1 .27e-04*ND  1  *LA2; 

KRDOT=-3.37e-05*NDl*LA3; 

KP=-1.10e-02*NDl*LA2; 

KPOMEGA=0; 

KP=KP+KPOMEGA*OMEGA; 

KV=3.055e-03*NDl*L; 

KVOMEGA=0; 

KV=KV+KVOMEGA*OMEGA; 

KR=-8.41e-04*NDl*LA2; 

KROMEGA=0; 

KR=KR+KROMEGA*OMEGA; 

flag=0; 
forXG=0:0.01:2, 

flag=flag+l; 

xg(flag)=XG; 

a=KX-KPDOT;  b=KP*U;  e=KV*U; 
f=KRDOT;  i=YP*U;  j=M-YVDOT;  k=YV*U;      ' 
l=XG*M-YRDOT;  m=U*(YR-M);  o=NPDOT; 
p=NP*U;  q=-XG*W;  r=XG*M-NVDOT; 
w=U*(NR-XG*M);  x=NV*U;  u=IZZ-NRDOT; 

al=-u*MA2; 

a2=-u*M*YPDOT-u*M*KVDOT+M*o*l+f*r*M; 

a3=a*j*u-a*l*r-u*KVDOT*YPDOT+KVDOT*o*l+f*r*YPDOT-f*o*j; 

b  1  =(w*  (MA2))+(r*U*  (MA2)); 

b2=-M*e*u-M*u*i+w*M*YPDOT+w*KVDOT*M-o*m*M+p*l*M- 

f*x*M+r*M*U*YPDOT-o*j*M*U+r*KR*U*M; 

b3=-e*u*YPDOT+e*o*l-u*i*KVDOT+w*KVDOT*YPDOT- 

KVDOT*o*m+KVDOT*p*l-a*j*w-a*k*u+a*l*x+a*r*m-b*j*u+b*l*r+f*r*i- 

f*x*YPDOT+o*k*f-f*p*j+r*KR*U*YPDOT; 

b2=b2+M*OMEGA*j*u-M*OMEGA*l*r; 

cl=-x*(MA2)*U; 

c2=r*i*M*U-x*M*U*YPDOT-x*KR*U*M+o*k*M*U-p*j*M*U+j*u*W-l*r*W- 

q*l*M+w*i*M-m*p*M+e*w*M; 

c3=a*k*w-a*m*x+b*j*w+b*k*u-b*l*x-b*r*m+r*i*KR*U- 

x*KR*U*YPDOT+o*k*KR*U-p*j*KR*U-q*l*KVDOT+w*i*KVDOT-m*p*KVDOT- 

e*u*i+e*w*YPDOT-e*o*m+e*p*l+f*q*j-f*x*i+f*p*k; 
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c2=c2-M*OMEGA*j*W-M*OMEGA*k*u-M*OMEGA*l*x+M*OMEGA*r*m; 

d2=q*j*M*U-x*i*M*U+p*k*M*U+q*m*M-j*w*W-k*u*W+l*x*W+r*m*W; 

d3=q*j*KJl*U-x*i*K^*U+p*k*KR*U-f*q*k+q*m*KVDOT-e*q*l+e*w*i-e*m*p- 

b*k*w+b*m*x; 

d2=d2+M*OMEGA*k*w-M*OMEGA*U*(KR-M)*x; 

e2=k*w*W-m*x*W-q*k*M*U; 
e3=e*q*m-q*k*KR*U; 

f5=(-e2*(blA2))+bl*cl*d2; 

f4=((b2*c  1  +b  1  *c2)*d2)+b  1  *c  1  *d3-al  *(d2A2)-e3  *(b  1  A2)-2*e2*b  1  *b2; 

f3=-a2*(d2A2)-2*d3*d2*al-2*e3*bl*b2- 

e2*((b2A2)+2*bl*b3)+d2*(b3*cl+b2*c2+bl*c3)+d3*(b2*cl+bl*c2); 

f2=-e3*((b2A2)+2*bl*b3)-2*e2*b2*b3+d2*(b3*c2+b2*c3)+d3*(b3*cl+b2*c2+bl*c3> 

a3*(d2A2)-2*d3*d2*a2-al  *(d3A2); 

fl=b3*c3*d2+d3*(b3*c2+b2*c3)-2*e3*b2*b3-e2*(b3A2)-2*d3*d2*a3-a2*(d3A2); 

f0=b3*c3*d3-a3*(d3A2)-e3*(b3A2); 

coef=[f5f4f3f2flf0]; 
ZG=roots(coef); 

tot(flag)=ZG(5,l); 
end 
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C      PROGRAM  HOPF_lNEW.FOR 

C      EVALUATION  OF  HOPF  BIFURCATION  FORMULAS  USING  THE  SUBOFF 

C      SUBMARINE  MODEL 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

REAL*8  L,IYY,M,YPDOT,YVDOT,NDl,YRDOT 

REAL*8YP,YV,YR,NPDOT,NVDOT,NRDOT,NP,NV,NR 

REAL*  8  GAMA,U,KVDOT,KRDOT,KPDOT,D 

REAL* 8  E0,E  1  ,E2,E3 ,E4,XG,ZG,KR,KP,KV,XG  1  ,ZG  1 
C 

REAL*8  Ml  1,M12,M13,M14,M21,M22,M23 

REAL*8M24,M31,M32,M33,M34,M41,M42,M43,M44 

REAL*8  Nl  1,N12,N13,N14,N21,N22,N23,N24 

REAL*8N31,N32,N33,N34,N41,N42,N43,N44 

REAL*8L11,L12,L13,L14,L21,L22,L23,L24,L31 

REAL*8L32,L33,L34,L15,L16,L17,L25,L26,L27,L35 

REAL*8L36,L37,L1A,L2A,L3A,L4A,L5A,L6A,L7A,L8A,L9A 

REAL*8  L10AJL1 1A,L12A,L1,L2,L3,L4,L5,L6,L7 

REAL*8  L8,L9,R1 1,R12,R13,R14,R21,R22,R23,R24 

REAL*8  P11,P12,P13,P21,P22,P23 

DOUBLE  PRECISION  THETA 
C 

DIMENSION  F(4,4),T(4,4),TINV(4,4),FV  1  (4),IV  1  (4), YYY(4,4) 

DIMENSION  WR(4),WI(4),TSAVE(4,4),TLUD(4,4),IVLUD(4),SVLUD(4) 

DIMENSION  AS AVE(4,4), A2(4,4),XL(  1 8),HT(  1 8),ZGG(202),FF(4,4) 

DIMENSION  VEC0(  1 8), VEC 1  ( 1 8), VEC2(  1 8), VEC3(  1 8), VEC4(  1 8),XGG(202) 

INTEGER  I,J,K 


C 
C 


OPEN  (20,FILE='HOPF25.RES') 
.     OPEN  (2 1  ,FILE='DAT25.DAT',STATUS=,OLD') 
C       OPEN  (23,FILE='HOPFl.RES,,STATUS='OLD') 
C 

WEIGHT=  12000.0 

IXX    =1760.0 

IYY    =9450.0 

IZZ    =10700.0 

EXZ    =0.0 

IXY    =0.0 

IYZ    =0.0 

L      =17.425 

RHO    =1.94 

CD     =0.5 
C      CD     =0.5*CD*RHO    ??? 

G      =32.2 
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XB     =0.0 

ZB     =0.0 

YG     =0.0 

YB     =0.0 

YDELTAR=0.0 

DELTAR=0.0 

NDELTAR=0.0 

NPROP=0.0 

M  =WEIGHT/G 

B    =  WEIGHT 

W=WEIGHT 

WRITE  (*,*)  '  ENTER  U" 
READ  (*,*)  U 

WRITE  (*,*) '  ENTER  THETA  (DEGREES)' 
READ  (*,*)  THETA 
C       U    =5.0 

ND1    =  0.5*RHO*L**2 
C       THETA=5 

THETA=THETA*PI/1 80 

OMEGA=U*DTAN(THETA) 
C 
C      DEFINE  HYDRODYNAMIC  COEFFICIENTS 

YPDOT=1.270E-04*ND1*L**2 

YVDOT=-5.550E-02*ND1*L 

YRDOT=1.240E-03*ND1*L**2 

YP=3.055E-03*ND1*L 

YPOMEGA=0.D0 

YP=YP+YPOMEGA*OMEGA+M*OMEGA 

YV=-9.310E-02*ND1 

YVOMEGA=0.D0 

YV=YV+YVOMEGA*OMEGA 

YR=-5.940E-02*ND1*L 

YROMEGA=0.D0 

YR=YR+YROMEGA*OMEGA 
C 

NPDOT=-3.370E-05*ND1*L**3 

NVDOT=1.240E-03*ND1*L**2 

NRDOT=-3 .400E-03  *ND  1  *L*  *  3 

NP=-8.405E-04*ND1*L**2 

NPOMEGA=0.D0 

NP=NP+NPOMEGA*OMEGA 

NV=- 1 .484E-02*ND  1  *L 

NVOMEGA=0.D0 

NV=NV+NVOMEGA*OMEGA 

NR=-1.640E-02*ND1*L**2 
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NROMEGA=0.D0 
NR=NR+NROMEGA*OMEGA 


C 
C 
C 


KPDOT=-1.01E-03*ND1*L**3 

KVDOT=1.27E-04*ND1*L**2 

KRDOT=-3.37E-05*ND1*L**3 

KP=-1.10E-02*ND1*L**2 

KPOMEGA=0.D0 

KP=KP+KPOMEGA*OMEGA 

KV=3.055E-03*ND1*L 

KVOMEGA=0.D0 

KV=KV+KVOMEGA*OMEGA 

KR=-8.41E-04*ND1*L**2 

KROMEGA=0.D0 

KR=KR+KROMEGA*OMEGA 

DEFINE  THE  LENGTH  X  AND  HEIGHT  H  TERMS  FOR  THE  INTEGRATION 
ALL  IN  FEET. 


XL(  1 

XL(2 

XL(3 

XL(4 

XL(5 

XL(6 

XL(7 

XL(8 

XL(9 

XL(10 

XL(11 

XL(12 

XL(13 

XL(14 

XL(15 

XL(16 

XL(17 

XL(18 

HT(  1 
HT(2 
HT(3 
HT(4 
HT(5 
HT(6 
HT(7 
HT(8 


=-105.9/12.0 
=-104.3/12.0 
=-99.3/12.0 
=-94.3/12.0 
=-87.3/12.0 
=-76.8/12.0 
=-66.3/12.0 
=-55.8/12.0 
=72.7/12.0 
=79.2/12.0 
=83.2/12.0 
=87.2/12.0 
=91,2/12.0 
=95.2/12.0 
=99.2/12.0 
=101.2/12.0 
=102.1/12.0 
=103.2/12.0 

=  0.000 
=  2.28/12.0 
=  8.24/12.0 
=  13.96/12.0 
=  19.76/12.0 
=  25.1/12.0 
=  29.36/12.0 
=  31.85/12.0 
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HT(  9)=  31.85/12.0 
HT(10)=  30.00/12.0 
HT(11)=  27.84/12.0 
HT(12)=  25.12/12.0 
HT(13)=  21.44/12.0 
HT(14)=  17.12/12.0 
HT(15)=  12.0/12.0 
HT(16)=  9.12/12 
HT(  17)=  6.72/12 
HT(18)=0.00 

DO  104  K=  1,18 
VEC0(K)=HT(K) 
VEC 1  (K)=XL(K)*HT(K) 
VEC2(K)=XL(K)*XL(K)*HT(K) 
VEC3(K)=XL(K)*XL(K)*XL(K)*HT(K) 
VEC4(K)=XL(K)*XL(K)*XL(K)*XL(K)*HT(K) 
104  CONTINUE 

CALL  TRAP(18,VEC0,XL,E0) 

CALL  TRAP(  1 8,VEC  1  ,XL,E  1 ) 

CALL  TRAP(  1 8,  VEC2,XL,E2) 

CALL  TRAP(  1 8,VEC3,XL,E3) 

CALL  TRAP(  1 8,VEC4,XL,E4) 

GAMA=0.001 


C       READ  THE  CRITICAL  VALUES  FOR  XG  AND  ZG  FROM  FILE  DATA  1  .DAT 
C 

XGG(1)=0.0 

ZGG(1)=0.0 16358083 
DO  1  1=2,202 

READ  (21,*)XG,ZG 
C  WRTTE(*,*)XG,ZG 

XGG(I)=XG 

ZGG(I)=ZG 


C-- 


C        DETERMINE  [F]  COEFFICIENTS 

C 

D=((M-YVDOT)*(IZZ-NRDOT)*(IXX-KPDOT)) 

&-((M-YVDOT)*(IXZ-KRDOT)*(-IXZ-NPDOT» 

&-((M*XG-NVDOT)*(M*XG-YRDOT)*(IXX-KPDOT)) 

&+((M*XG-NVDOT)*(IXZ-KRDOT)*(-M*ZG-YPDOT)) 

&+((-M*ZG-KVDOT)*(M*XG-YRDOT)*(-DCZ-NPDOT)) 
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&-((-M*ZG-KVDOT)*(IZZ-NRDOT)*(-M*ZG-YPDOT)) 

All=((IZZ-NRDOT)*(IXX-KPDOT))-((IXZ-KRDOT)*(-IXZ-NPDOT)) 

A12=((-M*XG+YRDOT)*(IXX-KPDOT))+((IXZ-KPvDOT)*(-M*ZG-YPDOT)) 

A13=((M*XG-YRDOT)*(-IXZ-NPDOT))-((IZZ-NRDOT)*(-M*ZG-YPDOT)) 

A21=((-M*XG+NVDOT)*(IXX-KPDOT))+((-M*ZG-KVDOT)*(-IXZ-NPDOT)) 

A22=((M-YVDOT)*(IXX-KPDOT))-((-M*ZG-KVDOT)*(-M*ZG-YPDOT)) 

A23=((-M+YVDOT)*(-IXZ-NPDOT))+((M*XG-NVDOT)*(-M*ZG-YPDOT)) 

A31=((M*XG-NVDOT)*(IXZ-KRDOT))-((-M*ZG-KVDOT)*(IZZ-NRDOT)) 

A32=((-M+YVDOT)*(IXZ-KPxDOT))+((-M*ZG-KVDOT)*(M*XG-YRDOT)) 

A33=((M-YVDOT)*(IZZ-NRDOT))-((M*XG-NVDOT)*(M*XG-YRDOT)) 


C-- 


C        EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 
C 

F(1,1)=(A11*YV*U+A12*NV*U+A13*KV*U)/D 

F(1,2)=(A11*(YR*U-M*U)+A12*(-M*XG*U+NR*U)+A13*(M*ZG*U+KR*U))/D 
F(  1 ,3)=(A  1 1  *  YP*U+A  1 2*NP*U+A  1 3  *(KP*U-M*ZG*OMEGA))/D 
F(l,4)=(All*(W-B)*DCOS(THETA)+A12*XG*(W-XB)*B*DCOS(THETA)+ 
&A13*(-ZG*W+ZB)*B*DCOS(THETA))/D 
F(2, 1  )=(A2 1  *  YV*U+A22*NV*U+A23*KV*U)/D 

F(2,2)=(A21*(YR*U-M*U)+A22*(-M*XG*U+NR*U)+A23*(M*ZG*U+KR*U))/D 
F(2,3)=(A2 1  *  YP*U+A22*NP*U+A23*(KP*U-M*ZG*OMEGA))/D 
F(2,4)=(A21*(W-B)*DCOS(THETA)+A22*(XG*W-XB*B)*DCOS(THETA)+ 
&A23*(-ZG*W+ZB*B)*DCOS(THETA))/D 
F(3,1)=(A31*YV*U+A32*NV*U+A33*KV*U)/D 

F(3,2)=(A3 1  *(YR*U-M*U)+A32*(-M*XG*U+NR*U)+A33*(M*ZG*U+KR*U))/D 
F(3,3)=(A3 1  *  YP*U+A32*NP*U+A33*(KP*U-M*ZG*OMEGA))/D 
F(3,4)=(A3 1  *(W-B)*DCOS(THETA)+A32*(XG*W-XB*B)*DCOS(THETA)+ 
&A33*(-ZG*W+ZB*B)*DCOS(THETA))/D 
F(4,l)=0.0 
F(4,2)=0.0 
F(4,3)=l-0 
F(4,4)=0.0 
C 

DOHK=l,4 
D0  12J=1,4 
ASAVE(K,J)=F(K,J) 
12      CONTINUE 
11    CONTINUE 

CALL  RG(4,4,F,  WR.WL  1 ,  YYY,IV  1  ,FV  1  ,IERR) 
WRITE(23,1007)WR(1),WR(2),WR(3),WR(4) 
CALL  DSOMEG(IEV,WR,WI,OMEGA,CHECK) 
C        WRITE  (*,*)  CHECK 
C         WRITE  (*,*)  WR(2) 
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C  WRITE  (*,*)  WI(4) 
OMEGA0=OMEGA 
D0  5J=1,4 

T(J,l)=YYY(J,ffiV) 

T(J,2)=-YYY(J,ffiV+l) 

5  CONTINUE 
IF(IEV.EQ.1.0)GOTO  13 
IF  (IEV.EQ.2.0)  GO  TO  18 
IF  (IEV.EQ.3.0)  GO  TO  14 

C         STOP  3004 
14    D0  6J=1,4 

T(J,3)=YYY(J,1) 
T(J,4)=YYY(J,2) 

6  CONTINUE 
GOTO  17 

18  D0  19J=1,4 
T(J,3)=YYY(J,1) 
T(J,4)=YYY(J,4) 

19  CONTINUE 
GOTO  17 

13    DO  16  J=l,4 

T(J,3)=YYY(J,3) 
T(J,4)=YYY(J,4) 

16  CONTINUE 

17  CONTINUE 
C 

C        NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 
C 

CALL  NORMAL(T) 
C 

C        INVERT  TRANSFORMATION  MATRIX 
C 

D0  2K=1,4 
D0  3J=1,4 
TINV(K,J)=0.0 
TSAVE(K,J)=T(K,J) 
3      CONTINUE 
2    CONTINUE 

C  ALL  DLUD(4,4,TS  AVE,4,TLUD,IVLUD) 
C        D0  4J=1,4 

C  IF  (IVLUD(J).EQ.O)  STOP  3003 

C     4    CONTINUE 

CALL  DILU(4,4,TLUD,IVLUD,SVLUD) 
D0  8K=1,4 
D0  9J=1,4 
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TINV(K,J)=TLUD(K,J) 
9   CONTINUE 
8  CONTINUE 
C 

C   CHECK  Inv(T)*A*T 
C 

CALL  MULT(TINV,ASAVE,T,A2) 
C 

P1=A2(3,3) 

P2=A2(4,4) 
PEIG1=P1 
PEIG2=P2 
C 

C        DEFINITION  OF  Nij 
C 

N11=TINV(1,1) 
N12=TINV(1,2) 

N13=TINV(1,3) 
N14=TINV(1,4) 

N21=TINV(2,1) 

N22=TINV(2,2) 
N23=TINV(2,3) 
N24=TINV(2,4) 
N31=TINV(3,1) 
N32=TINV(3,2) 

N33=TINV(3,3) 

N34=TINV(3,4) 

N41=TINV(4,1) 

N42=TINV(4,2) 

N43=TINV(4,3) 

N44=TINV(4,4) 
C 

C        DEFINITION  OF  Mij 
C 

M11=T(1,1) 

M12=T(1,2) 

M13=T(1,3) 
M14=T(1,4) 
M21=T(2,1) 

M22=T(2,2) 

M23=T(2,3) 
M24=T(2,4) 
M31=T(3,1) 
M32=T(3,2) 
M33=T(3,3) 
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M34=T(3,4) 

M41=T(4,1) 

M42=T(4,2) 

M43=T(4,3) 

M44=T(4,4) 
C 

C        DEFINITION  OF  Lij 
C 

L1=YG*((M31**2)+(M21**2)) 

L2=IXY*(M31**2)+IYZ*M31*M21+YG*M21*M11 

L3=IXY*M31*M21+IYZ*(M21**2)+M*YG*M11*M31+YG*W*(M41**2)- 
YB*B*(M41** 

&2) 

L4=YG*(2.0*M3 1  *M32+2.0*M2 1  *M22) 

L5=2.0*IXY*M3 1  *M32+IYZ*(M3 1  *M22+M32*M2 1  )+YG*(M2 1  *M  1 2+M22*M  1 1 ) 

L6=KY*(M3 1  *M22+M32*M2 1  )+2.0*IYZ*M2 1  *M22+M*  YG*(M  1 1  *M32+M  1 2*M3 1 
)+2. 

&0*(YG*W-YB*B)*M41*M42 

L7=YG*((M32**2)+(M22**2)) 

L8=KY*(M32**2)+IYZ*M32*M22+YG*M22*M12 

L9=rXY*M32*M22+IYZ*(M22**2)+M*YG*M12*M32+(YG*W- 
YB*B)*(M42**2) 
C 
C 

L15=(A11/D)*L1+(A12/D)*L2-(A13/D)*L3 

L16=(A11/D)*L4+(A12/D)*L5-(A13/D)*L6 

L17=(A11/D)*L7+(A12/D)*L8-(A13/D)*L9 

L25=(A21/D)*L1+(A22/D)*L2-(A23/D)*L3 

L26=(A21/D)*L4+(A22/D)*L5-(A23/D)*L6 

L27=(A21/D)*L7+(A22/D)*L8-(A23/D)*L9 

L35=(A31/D)*L1+(A32/D)*L2-(A33/D)*L3 

L36=(A31/D)*L4+(A32/D)*L5-(A33/D)*L6 

L37=(A31/D)*L7+(A32/D)*L8-(A33/D)*L9 

C 

C=CD/(6.0*GAMA) 

L'1A=C*(E0*(M11**3)+3.0*E1*(M11**2)*M21+3.0*E2*M11*(M21**2)+E3*(M21 
&**3))-((  1 .0/6.0)*(W-B)*DCOS(THETA))*(M4 1  **3) 

L2A=C*(E1*(M11**3)+3.0*E2*(M11**2)*M21+3.0*E3*M11*(M21**2)+E4*(M21 
«&**3))-((1.0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M41**3) 
L3A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M41**3) 
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L4A=C*(3.0*E0*(Mll**2)*M12+3.0*El*((Mll**2)*M22+2*Mll*M12*M21)+3.0 
&*E2*(M12*(M21**2)+2.0*M21*M22*Mll)+3.0*E3*(M21**2)*M22)-((1.0/6.0) 
&*(W-B)*DCOS(THETA))*3.0*(M41**2)*M42 

L5A=C*(3.0*El*(Mll**2)*M12+3.0*E2*((Mll**2)*M22+2*Mll*M12*M21)+3.0 
&*E3*(M12*(M21**2)+2*M21*M22*M11)+3.0*E4*(M21**2)*M22)-((1.0/6.0)*(X 
&G*W-XB*B)*DCOS(THETA))*3.0*(M41**2)*M42 
L6A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*(M41**2)*M42 

L7A=C*(3*E0*M11*(M12**2)+3.0*E1*((M12**2)*M21+2*M11*M12*M22)+3.0*E 
&2*((M22**2)*M11+2.0*M21*M22*M12)+3.0*E3*M21*(M22**2))-((1.0/6.0)*(W 
&-B)*DCOS(THETA))*3.0*M41*(M42**2) 

L8A=C*(3.0*El*Mll*(Mi2**2)+3.0*E2*((M12**2)*M21+2.0*Mll*M12*M22)+3 
&.0*E3*((M22**2)*M11+2*M21*M22*M12)+3.0*E4*M21*(M22**2)) 
L9A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*M41*(M42**2) 
L10A=-C*(E0*(M12**3)+3.0*E1*(M11**2)*M21+3.0*E2*M12*(M22**2)+E3*(M 
&22**3))-((1.0/6.0)*(W-B)*DCOS(THETA))*(M42**3) 

L11A=-C*(E1*(M12**3)+3.0*E2*(M11**2)*M21+3.0*E3*M12*(M22**2)+E4*(M 
&22**3))-((l-0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M42**3) 
L12A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M42**3) 


Lll  = 
LI  2= 
L13= 
L14= 

L21  = 
L22= 
L23= 
L24= 

L31  = 
L32= 
L33= 
L34= 

Rll= 
R12= 
R13= 
R14= 
R21= 
R22= 
R23= 


-A11/D)*L1A+(A12/D)*L2A+(A13/D)*L3A 
-Al  1/D)*L4A+(A12/D)*L5A+(A13/D)*L6A 
-A11/D)*L7A+(A12/D)*L8A+(A13/D)*L9A 
-Al  1/D)*L10A+(A12/D)*L1 1A+(A13/D)*L12A 

-A21/D)*L1A+(A22/D)*L2A+(A23/D)*L3A 
-A21/D)*L4A+(A22/D)*L5A+(A23/D)*L6A 
-A21/D)*L7A+(A22/D)*L8A+(A23/D)*L9A 
-A21/D)*L10A+(A22/D)*L11A+(A23/D)*L12A 

-A31/D)*L1A+(A32/D)*L2A+(A33/D)*L3A 
-A3 1/D)*L4A+(A32/D)*L5A+(A33/D)*L6A 
-A31/D)*L7A+(A32/D)*L8A+(A33/D)*L9A 
-A31/D)*L10A+(A32/D)*L11A+(A33/D)*L12A 

N11*L11)+(N12*L21)+(N13*L31) 
N11*L12)+(N12*L22)+(N13*L32) 
Nl  1*L13)+(N12*L23)+(N13*L33) 
Nl  1*L14)+(N12*L24)+(N13*L34) 
N21*L11)+(N22*L21)+(N23*L31) 
N21*L12)+(N22*L22)+(N23*L32) 
N21*L13)+(N22*L23)+(N23*L33) 
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R24=(N21*L14)+(N22*L24)+(N23*L34) 
C 

P11=(N11*L15)+(N12*L25)+(N13*L35) 
P12=(N11*L16)+(N12*L26)+(N13*L36) 
P13=(N11*L17)+(N12*L27)+(N13*L37) 
P21=(N21*L15)+(N22*L25)+(N23*L35) 

P22=(N21*L16)+(N22*L26)+(N23*L36) 
P23=(N21*L17)+(N22*L27)+(N23*L37) 
C 

C        EVALUATE  DALPHA  AND  DOMEGA 
C 

ZGR  =ZGG(I) 

ZGL=ZGG(M) 

ZG1  =ZGR 

XG1  =XGG(I) 
C 

CALL  FMATRIX(ZG1,XG1,FF,U,THETA) 
C 

CALL  RG(4,4,FF,  WR,  WI,0,  YYY,IV  1  ,FV  1 JERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 
C 

ZG1  =ZGL 

XG1  =XGG(I-1) 
C 

CALL  FMATRIX(ZG1,XG1,FF,U,THETA) 
C 

CALL  RG(4,4,FF,WR,WI,0,  YYY,IV  1  ,FV  1  ,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 
C 

DALPHA=(ALPHR-ALPHL)/(ZGR-ZGL) 

DOMEGA=(OMEGR-OMEGL)/(ZGR-ZGL) 
C 

C        EVALUATION  OF  HOPF  BIFURCATION  COEFFICIENTS 
C 

COEF 1  =(  1 .0/8 .0)  *  (3 .0*R  1 1 +R 1 3+R22+3 .0*R24) 
COEF2=(1.0/8.0)*(3.0*R11+R23-R12-3.0*R14) 

C  AMU2  =-COEFl/(8.0*DALPHA)     ???????? 

C  BETA2=0.25*COEF1  ??????? 

C  TAU2  =-(COEF2-DOMEGA*COEFl/DALPHA)/(8.0*OMEGAO) 

C  PER  =2.0*3.1 41 5927/OMEGA0 
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C       PER  =PER*U/L 

WRITE  (20,200 1)XG,ZG,COEF1,DALPHA,OMEGAO,PEIG1,PEIG2 
C       WRITE  (20,200 1  )COEF  1 

1  CONTINUE 
C 

STOP 

2001  FORMAT  (7E14.5) 
4001  FORMAT  (3F15.5) 
1007  FORMAT  (4E  14.5) 
END 
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